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Abstract 

We describe a uniformly fast algorithm for generating points x 
uniformly in a hypercube with the restriction that the difference be- 
tween each pair of coordinates is bounded. We discuss the quality of 
the algorithm in the sense of its usage of pseudo-random source num- 
bers, and present an interesting result on the correlation between the 
coordinates. 
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1 Introduction 



In this paper we shall discuss the problem of generating sets of points x = 
(xi, X2, ■ ■ • , x m ) inside an m-dimensional hypercube with an additional re- 
striction. The points x are required to satisfy the conditions 

\xk\ < 1 , \xk — xi\ < 1 for all k, I . (1) 

These conditions define a m-dimensional convex polytope P. The reason for 
tackling this problem is the following. In a recently developed Monte Carlo 
algorithm, SARGE [jl], we address the problem of generating configurations of 
four-momenta pf, % = 1, 2, . . . , n of n massless partons at high energy, with 
a distribution that has, as much as possible, the form of a so-called QCD 
antenna: 

Ski = (Pk+Pi) 2 , 



Sl2S23 s 34 ' ' ' S n -i >n S n i 

where is the invariant mass squared of partons k and I, with the addi- 
tional requirement that the total invariant mass squared of all the partons is 
fixed to s, and every s^i (also those not occurring explicitly in the antenna) 
exceeds some lower bound so'- in this way the singularities of the QCD matrix 
elements are avoided. The SARGE algorithm has a structure that is, in part, 
similar to the RAMBO algorithm ||, where generated momenta are scaled so 
as to attain the correct overall invariant mass. Obviously, in SARGE this is 
more problematic because of the Sq cut, but one should like to implement 
this cut as far as possible. Note that out of the n(n — l)/2 different s^i, n 
occur in the antenna, and each of these must of course be bounded by s 
from below and some % < s from above. The scale- invariant ratios of two 
of these masses are therefore bounded by 

il < S JL < s Ji , ( 2 ) 

% Ski Sq 

The structure of the SARGE algorithm is such ]I| that there are m = 2n — 4 
of these ratios to be generated. By going over to variables 

£(...) =log(s ij /s k i)/log( y S M /so) , 

and inspecting all ratios that can be formed from the chosen m ones, we 
arrive at the condition of Eq.(|l|). Note that, inside SARGE, a lot of internal 
rejection is going on, and events satisfying Eq.(|l|) may still be discarded: 
however, if Eq. Ql|) is not satisfied, the event is certainly discarded, and it 
therefore pays to include this condition from the start. 
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2 The algorithm 



The most straightforward way of implementing is of course the following: 
generate Xk, k = 1, . . . , to by Xk <— 2p — 1, and reject if the conditions 
are not met. Here and in the following, each occurrence of p stands for a 
call to a source of iid uniform pseudo-random numbers between in [0, 1). The 
drawback of this approach is that the efficiency, i.e. the probability of success 
per try, is given by 2~ m V m (P) (where V m (P) is the volume of the polytope 
P) and becomes very small for large to, as we shall see. 

To compute the volume V m (P) we first realize that the condition \xk — 
x\\ < 1 is only relevant when Xk and x\ have opposite sign. Therefore, we 
can divide the x variables in to — k positive and k negative ones, so that 
i 

V m ,k(P) = J dy 1 dy 2 ■ ■ ■ dy k dxk+idx k+ 2 • • • dx m 9 (l - rnaxi, - maxj/j 
o 

v^p) - tj^^niP) , (?) 

where we have written y k = —Xk- By symmetry we can always relabel the 
indices such that x m = max^Xj and y\ = maxj-y^. The integrals over the 
other x's and y's can then easily be done, and we find 

1 l-2/l 

V mik (P) = k(m -k)j dywt 1 J dx m xZ' k - 1 

o o 

= k]dy iy \-\l-yr- k = kKm ~ k) - , (4) 
o 

and hence 

V m {P)=m + l . (5) 

The efficiency of the straightforward algorithm is therefore equal to (to + 
l)/2 m , which is less than 3% for n larger than 6. 

We have given the above derivation explicitly since it allows us, by work- 
ing backwards, to find a rejection-free algorithm with unit efficiency. The 
algorithm is as follows: 

1. Choose a value for k. Since each k is exactly equally probably we simply 
have 

k <- [{m + l)p\ . 
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2. For k = we can simply put 

Xi <- p ,i = 1, . . . ,m , 

while for k = m we put 

Xi < p ,i = 1, ...,ra . 

3. For < k < m, y± has the unnormalized density — y\) m ~ k 
between and 1. An efficient algorithm to do this is Cheng's rejection 
algorithm BA for beta random variates (cf. 0)[], but the following also 
works: 

/ k \ /m-k+l \ ^ 

vi< log [YI p) , v 2 < log n p , yi 



\i=l 



j=i I Vi+V 2 



The variable x m has unnormalized density x™ between and 1 — yi 
so that it is generated by 

x m ^(l-2/i)p 1/(m ~ fc) • 
The other x's are now trivial: 



x%< yi , Xi <— xtp, £ = 2,3, , 

Xi <— x m p, i — k + 1, k + 2, . . . , m — 1 . 

Finally, perform a random permutation of the whole set (xi, x 2 , ■ ■ ■ , x m ). 



3 Computational complexity 

The number usage S, that is, the expected number of calls to the random 
number source p per event can be derived easily. In the first place, 1 number 
is used to get k for every event. In a fraction 2/(m + 1) of the cases, only m 
calls are made. In the remaining cases, there are k + (m — k + 1) = m + l 
calls to get yi, and 1 call for all the other x values. Finally, the simplest 

1 There is an error on page 438 of [fj], where "V <— A _1 J7i(l — Ui)^ 1 " should be replaced 
by U V<- X-Hogpiil-Ux)- 1 }". 
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permutation algorithm calls m — 1 times [[|. The expected number of calls 
is therefore 

2m m— 1. „ , „ , , 3m 2 — m + 2 . , 

5 = 1 + - + m + l + (m-1 + m-1 = . 6 

m + 1 m + l v v 7 v " m+1 w 

For large m this comes to about 3m — 1 calls per event. Using a more 
sophisticated permutation algorithm would use at least 1 call, giving 

2m* in — 1 

S = l + + Am + 1 + (m - 1) + (1)) = 2m . (7) 

m + 1 m+ l v v ) \ n \ i 

We observed that Cheng's rejection algorithm to obtain y\ uses about 2 calls 
per event. Denoting this number by C the expected number of calls becomes 

am' + tc-pm-c + a 

m+1 v 7 

for the simple permutation algorithm, while the more sophisticated one would 
yield 

m' + tC + qm-C+l 

m+1 w 
We see that in all these cases the algorithm is uniformly efficient in the sense 
that the needed number of calls is simply proportional to the problem's 
complexity m, as m becomes large. An ideal algorithm would of course 
still need m calls, while the straightforward rejection algorithm rather has 
S = m2 m / (m + 1) ~ 2 m expected calls per event. 

In the testing of algorithms such as this one, it is useful to study expec- 
tation values of, and correlations between, the various X{. Inserting either X{ 
or XiXj in the integral expression for V(P), we found after some algebra the 
following expectation values: 

Efe ) = ° • E «> = 6|S7TT) • ^=ufrTT) (10) 

so that the correlation coefficient between two different x's is precisely 1/2 
in all dimensions! This somewhat surprising fact allows for a simple but 
powerful check on the correctness of the algorithm's implementation. 

As an extra illustration of the efficiency, we present in the tables below the 
cpu-time (t cpu ) needed to generate 1000 points in an m-dimensional polytope, 
both with the algorithm presented in this paper (OURALG) and the rejection 
method (REJECT). In the latter, we just 
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1. put Xi <— 1p — 1 for i — 1, . . . , m; 

2. reject x if \xi — Xj \ > 1 for i = 1, . . . , m — 1 and j — % + 1, . . . , m. 

The computations were done using a single 333-MHz UltraSPARC-Hi pro- 
cessor, and the random number generator used was RANLUX on level 3. 





^cpu(scc) 


m 


OURALG 


REJECT 


2 


0.03 


0.01 


3 


0.03 


0.02 


4 


0.03 


0.04 


5 


0.04 


0.08 


6 


0.05 


0.17 


7 


0.06 


0.32 


8 


0.07 


0.67 


9 


0.08 


1.33 


10 


0.09 


2.76 



m 


OURALG 


REJECT 


11 


0.09 


5.15 


12 


0.10 


10.94 


13 


0.11 


21.71 


14 


0.12 


44.06 


15 


0.13 


87.90 


16 


0.14 


169.65 


17 


0.15 


336.67 


18 


0.16 


671.46 


19 


0.17 


1383.33 


20 


0.18 


2744.82 



For m = 2 and m = 3, the rejection method is quicker, but from m = 4 on, 
the cpu-time clearly grows linearly for the method presented in this paper, 
and exponentialy for the rejection method. 



4 Extension 

Let us, finally, comment on one possible extension of this algorithm. Suppose 
that the points x are distributed on the polytope P, but with an additional 
(unnormalized) density given by 

F(f) = flcos(^) , (11) 

so that the density is suppressed near the edges. It is then still possible to 
compute V mt k(P) for this new density: 

i 1-2/1 
V k , m (P) = k(m - k) j dyi cos (-y-) f dx m cos^' 
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yi \ k ~ x / x m \ m-k-l 

{1TX 

ax cos — 



/^ cos (t 

vo 



a sin sin 



fc-i 



2 m ~ 1 k 



r'" 



2 / V V 2 

| ds s fc/2-l(l_ a )(™-*)/2 



7f 



2\ m T(l + fc/2)r(l + (m - fc)/2) 



7T 



r(l + m/2) 



(12) 



where we used s = (sin (^f 1 )) • Therefore, a uniformly efficient algorithm 
can be constructed in this case as well, along the following lines. Using the 
Vfe >m , the relative weights for each k can be determined. Then s is generated 
as a (3 distribution. The generation of the other x's involves only manipula- 
tions with sine and arcsine functions. Note that, for large m, the weighted 
volume of the polytope P is 

v(p) = fry Mg): ml 

1 ' fr VW (f)! k\(m-k)\ 



so that a straightforward rejection algorithm would have number usage 

r x 

and a correspondingly decreasing efficiency 
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